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We generalize the imaginary chemical potential quantum Monte Carlo (QMC) method proposed 
by Dagotto et al. [Phys. Rev. B 41, R811 (1990)] to systems without particle-hole symmetry. The 
generalized method is tested by comparing the results of the QMC simulations and exact diagonal- 
ization on small Hubbard molecules, such as tetrahedron and truncated tetrahedron. Results of the 
application of the method to the Ceo Hubbard molecule are discussed. 
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I. INTRODUCTION 

Knowledge of the evolution of energy levels with dop- 
ing in strongly correlated systems is of significant impor- 
tance for understanding the physical mechanisms leading 
to their unconventional properties. For example, infor- 
mation about the changes in the ground state energy of 
model electron systems (such as Hubbard model) upon 
electron or hole doping may be used to confirm or dis- 
prove hypotheses about the origin of the pairing mecha- 
nism, eventually leading to superconductivity. As access- 
ing this information analytically usually requires the use 
of various approximations, numerical techniques are of- 
ten the only tools which can provide unbiased estimates 
for the observables of interest. 

We have recently applied the auxiliary field quantum 
Monte Carlo (AFQMC) on a Cm molecule to extract 
the electronic binding energies AFQMC has been 
widely used in Hubbard Hamiltonian simulations since 
its introduction by Blankenbecler et al. 0, Q > ^'^d its 
further development by Hirsch I] and White et al. j^. 
Being a finite-temperature technique, AFQMC does not 
allow easy access to the physical observables, not repre- 
sented by thermodynamic averages, such as energy gaps. 
A convenient procedure to extract this additional infor- 
mation from the AFQMC data was proposed by Dagotto 
et al. ,6J, who introduced imaginary chemical potentials 
in AFQMC simulations. It was then used to extract the 
charge gaps of the one-band Hubbard model on finite 
two-dimensional (2D) square lattices. 

In the present paper we generalize this formalism to 
systems without particle-hole symmetry, such as the 
tetrahedron, truncated tetrahedron and Ceo molecules. 
The canonical partition function ratios are obtained from 
the expansions of the AFQMC determinant ratios for a 
set of finite temperatures T, which are subsequently used 
to extract charge gaps at low temperatures. This general- 
ization results in the appearance of an extra phase factor 
in the expansion of the determinant ratios, which reduces 
to unity in the particle-hole symmetric systems. 

The rest of the paper is organized as follows. First, 
we briefly describe the imaginary chemical potential 
QMC (ICPQMC) formahsm. Then simulation results 
on some Hubbard molecules are presented to illustrate 



our method. The results are compared with the data ob- 
tained by exact diagonalization (ED) on small molecules 
and projector QMC (PQMC) on larger ones. 



II. METHODOLOGY 

We start with an expansion of the grand canonical par- 
tition function Z^ic{^^') in terms of canonical partition 
functions Zc{n) [g: 

N 

Zgc(a*) = Tre-'^(^-'^*) =6^^^^ J2 e^^"Zc(n), (1) 

where f3 — l/(fcsT) is the inverse temperature, /i is the 
chemical potential, n is the deviation of the particle num- 
ber from half-filling (positive or negative, denoting elec- 
tron or hole doping, respectively) in canonical ensemble, 
N is the electron number operator, and N is the num- 
ber of spatial lattice sites in the system. H is the usual 
one-band Hubbard Hamiltonian: 

iJ = - ^ tij{cl^Cja + h.c.) + n.^n^i - X! 

(2) 

The summation in the hopping term t.y is performed over 
all nearest-neighbor pairs of the Hubbard molecule. For 
the Ceo molecule we have set tij = t for the links between 
the pentagons and hexagons and tij = 1.2t for the links 
between hexagons. In all other cases tij was set equal to 
t for all links, with t used as an energy unit. U is the 
on-site Coulomb repulsion (Hubbard) term, and an extra 
diagonal term has been added to the Hamiltonian so that 
/i = corresponds to half-filling on bipartite lattices. 

Following Dagotto et al. 0, we analytically continue 
Eq. to imaginary chemical potential fj, iX, where 
A is real. Then the inverse Fourier transform of Eq. (jT]) 
gives 

fl /.27r//3 

Zc{n) = j- dX e-'^^("+^)ZGc(A* = *A). (3) 

Stt Jo 

In AFQMC the grand canonical partition function is 
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FIG. 1; Fit of real (squares) and imaginary (circles) parts 
of determinant ratios according to Eq. 1121 for a tetrahedron 
molecule. 



given by 

^gc(m) - II n det[l + BL{a)BL-i{a) ■ ■ ■ Bi{a)] 

{a} a = ±l 

= Y,detO{{a},^l)^detO{{a},^l)^, (4) 
{'^} 

where the fermion degrees of freedom have been traced 
out, and the Bi matrices are defined as 



Bi{a) 



-tij for ij nearest neighbors, 



'''^ 10 otherwise. 



(5) 
(6) 
(7) 



Here Ar is the imaginary time discretization interval and 
tanh^(7/2) = tanh(Art7/4). Thus, the original prob- 
lem of taking a trace over fermionic degrees of freedom 
has been replaced by a problem of tracing over auxiliary 
Ising variables cri{l), introduced at every space-time point 
{iJAr). Inserting Eq. into Eq. ©, we get 

r. ,27r//3 

^cW = / dAe-'^^^("+^)detO(M,iA)T 

xdctO(M,zA)4. (8) 

Dividing Eq. ^ by 2'gc(m = 0) from Eq. igj) yields 

Zcjn) ~ « /■^'^/'^ 

^GC(O) 



E^({'^}'0)f / dXe'^^'-A{pX) (9) 

{a} 



where 

A(/3A) = e 
and 



^ ^zl3\N det 0({o-}, det 0({cr}, zA)x 

detO(M,0)TdetO(M,0)x ' ^ ^ 



p. r . .^ ^ detO(M,0)TdetO({a},0)x 



o 

N 



O 



o 

N 



O 



-1.5 

-2 
-2.5 

-3 
-3.5 

-4 
-4.5 

-5 
-5.5 

-6 
-6.5 



3 
2.5 

2 
1.5 

1 

0.5 




T 1 1 r 

(a) Tetrahedron, U=2t 
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(b) Tetrahedron, U=2t 




FIG. 2: Fits of ICPQMC data for a tetrahedron molecule at 
low temperatures for electron (a) and hole (b) doping. NLF 
forms for (a) and (b) are given by Eq. 11141 and Eq. I15II . ED 
results are also shown for comparison. 



is the probability distribution for Zgc(O). Since there 
is always an energy gap above and below half-filling for 
any finite system, we expect the fermion determinants 
to be nearly A-independent at low temperatures. There- 
fore, we can generate the Ising field configurations {a} 
for A = and use these configurations to calculate system 
properties at A / 0. 

Similar to the expansion in Eq. we expect that the 
determinant ratio A(/3A) in Eq. can be expressed as 
complex Fourier series in the particle number n 

N 

A(/3A) = co(W) + E[c«(W) + c-„(W)]cos(/?An) 

AT 

+zE[cn({a}) -c_„({a})]sin(/3An), (12) 



where n {—n) represents a doping of n electrons (holes) 
with respect to half-filling. Eq. H12|l is real for sys- 
tems with particle-hole symmetry, since then c„({(t}) — 
c-n{{<j}). When we substitute Eq. 112|l back into Eq. 
©, we see that Zc{n) / Zqc{0) — {cn), where the aver- 
age (. . .) is over the Ising field configurations {a} gen- 
erated from the probability distribution P{{a},0). In 
the case of negative weight we replace P by its abso- 
lute value \P\ and include a sign 5 = P/\P\ in the av- 
erage: Zc{n)/ZGc{0) = {cnS)/{S). The average (. . .) 
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LF 


NLF 


ED (PQMC) 






Ai,o 


0.85(1) 


0.87(2) 


0.82843 


(U = 


2t) 


A_i,o 


-0.458(9) 


-0.358(9) 


-0.34949 






A2,l 


1.9(1) 


2.0(1) 


2.0 






A_2,-l 


-0.400(6) 


-0.397(6) 


-0.40466 


Cs 




A±i,o 


1.20(2) 


1.27(2) 


1.26224 


iu = 


4t) 


^^zb2 dzl 


1 44f3) 


1 5(2) 


1.27490 


Cl2 




Ai,o 


0.81(1) 


0.997(7) 


0.99596 


(C/ = 


2t) 


A-1,0 


0.041(1) 


0.10(1) 


0.07408 


Ceo 




Ai,o 


0.43(5) 


0.57(3) 


0.561(7)* 


([/ = 


4t) 


A_i,o 


0.9(1) 


0.88(4) 


0.86(2)* 



TABLE I: ICPQMC results on tetrahedron (C4), cube (Cs), 
truncated tetrahedron (C12) and Ceo molecules. A„.o = 
E{n)—E{0) is the energy difference between the ground states 
of the two fillings. Data marked by * are PQMC results, as 
described in Section fill Bl 



now refers to the probability distribution As we are 
interested in regions near half-filling, the sign problem 
does not limit the applicability of our method. 

At low temperatures the canonical partition function 
ratio will be dominated by A„^o = E{n) — E{Q), the 
energy difference between the ground states for the two 
fillings, and will take the form 



Zcjn) 
ZciO) 



(co) 



dn.oe 



as (3 



(13) 



where dn,o = dn/do, with (i„ being the degeneracy of the 
ground state at filling n. When there is an energy level 
close to the ground state (produced by the elementary 
excitations such as spin waves), we include it explicitly 
in the fitting expressions: 



Zcjn) 
ZciO) 



1 + /o,oe-^^^'w 
1 1 



as 



(14) 



Here /„^o — d^^/do, where d^^ is the degeneracy of the 
spin wave state at fiUing n, and A"^ = Esv,{n) — E{0) are 
the spin-wave gaps at filling n with respect to the ground 
state at half-filling. 

Based on the above discussion, we formulate the fol- 
lowing calculation procedure: 

1. Generate the Ising field configuration {a} accord- 
ing to the probability distribution P{{a}, 0) in the 
AFQMC simulation of ZociO)- 

2. Evaluate the average of the determinant ratio on 
the left-hand side of Eq. H12|) over the Ising field 
configurations for a set of A values. 

3. Fit the real and imaginary parts of Eq. (|12|l re- 
spectively to determine the average values of (c„), 
n = 0, ±1, ±2, . . .. From now on we will refer to 
these values simply as c„. 
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FIG. 3: Fits of the real (squares) and imaginary (circles) parts 
of determinant ratios according to Eq. (1121 for a Ceo molecule. 



4. Fit the canonical partition function ratios 
Zc(j>')/Zc{0) = Cn/cg, n = ±1,±2,... with low 
temperature canonical partition function ratio 
expressions ((T^ or lfn|) to obtain energy gaps. 
Below we will refer to these fits as linear (LF) and 
non- linear (NLF), respectively. 



III. APPLICATION 

ICPQMC simulations have been carried out on tetra- 
hedron, cube, truncated tetrahedron and Ceo molecules. 
For each molecule, we have run the simulations at tem- 
peratures f3t = 2.0, 2.5, • ■ • , 6.5, 7.0. The imaginary 
chemical potential iX was chosen so that < A < 7r//3, 
and we used a set of 20 evenly distributed A values in 
this range. We have numerically checked that the real 
part of Eq. H12() is an even function around A = 0, while 
the imaginary part is odd. Using this property, we have 
mapped out the determinant ratio data for A > 0. For 
the special particle-hole symmetric case, such as a 2D 
square lattice, we have tested our programs for the 2x2 
and 4x4 systems, reproducing the results of Ref. |^. 



A. Tetrahedron, Cube and Truncated Tetrahedron 

Fig- n displays the fit of the real and imaginary parts 
of Eq. H12|l for a tetrahedron molecule (C4). Similar fits 
were performed for data at all temperatures, generating 
a set of coefficients c„. Knowledge of these coefficients 
enables us to calculate the partition function ratios at 
various temperatures, and, at low temperatures, to lin- 
early fit the logarithm of these ratios to obtain the energy 
gaps and degeneracy ratios. Results of this procedure are 
presented in Fig. Inland Table We see that A2,i and 
A_2,-i from LF agree nicely with ED, while A_i_o does 
not. A likely cause of this discrepancy is the existence 
of a highly degenerate (d^^ = 9) spin wave energy level 
very close (A^^ = 0.14258i) to the ground state (do = 2) 
at half-filling. Results of the fits for electron (n = 1) and 



4 



hole (n = —1) doping are presented in Figs. |^ and[2l3 
respectively. The NLFs were obtained using form 1)14(1 
by fixing the spin wave and degeneracy parameters to 
the values found by ED: 



/i,o — 3, 
/o,o = 4.5, 



Al^ = 1.18268i, 
= 0.14258t. 



(15) 



In the case of hole doping the spin wave term propor- 
tional to /-i,o in the numerator of Eq. p4|l has been 
neglected. 

Examining Fig.|21we find that inclusion of spin waves 
does not make a substantial difference: LF and NLF 
curves are nearly overlapping [with a slight difference 
around pt = 2 for ln(Zc(l)/-^c(0))] and are very close 
to the exact values. This insensitivity to the fitting pa- 
rameters makes accurate extraction of energy gaps from 
partition function ratios difficult. 

Analogous simulations and LF/NLF fittings were per- 
formed for cube (Cs) and truncated tetrahedron (C12) 
molecules. The energy gaps measured using this proce- 
dure are summarized in Table The fitting procedure 
for Cgo is described below. 



B. 



Fig. O shows a fit of the real and imaginary parts of 
Eq. H12|l for a Ceo molecule. The imaginary parts of the 
determinant ratios are positive for A > 0, which is dif- 
ferent from the tetrahedron and truncated tetrahedron 
cases. The cause of this difference is the relative size 
of electron and hole gaps in the system. The larger the 
gaps, the smaller the corresponding canonical partition 
function ratios Zc{n) / Zc{0). For the tetrahedron and 
the truncated tetrahedron the electron gaps are larger 
than the hole gaps, so ci < c_i. In contrast, for the Cgo 
molecule the electron gap is smaller than the hole gap, 
so ci > c_i. Therefore, the relative size of the canoni- 
cal partition function ratios results in a positive imagi- 
nary part of the determinant ratios for A > in Fig. |3| 
due to Eq. (|12|) . Similar fits were done for other tem- 
peratures, and the resulting canonical partition function 
parameters c„ were obtained to calculate the canonical 
partition function ratios. Unfortunately, Eq. ((14|l con- 
tains too many fitting parameters to provide unique fits 
to the data. Therefore, we had to rely on the PQMC 
results (at U = At) for the gap values 



Ai,o 



0.561i, 
1.39i, 



ArJ 



1.54t, 



(16) 



and on the analysis of the molecular orbital energy level 
diagram (Fig. 3 of Ref. ilj) for the degeneracies and their 
ratios 



We were able to demonstrate the consistency of the 
results obtained by the two methods (ICPQMC and 
„ 1 
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(17) 



FIG. 4: Fits of ICPQMC data for a Ceo molecule at low 
temperatures for electron and hole doping, e. (h.) denotes 
electron (hole) doping. 



PQMC) by employing a following procedure. We first let 
all three gap values be free fitting parameters. Fitting 
results are then in agreement with the PQMC results, 
albeit with rather large uncertainties. Then we fill in the 
two PQMC gap values A^^ and [Eq. ITHJl], and let 
A-ti be the only free parameter. This procedure, in 
general, yields the gap values consistent with the ones 
previously obtained using PQMC method Q. Detailed 
results are presented in Table with representative fits 
shown in Fig. 21 



IV. CONCLUSION 

We have generalized the particle-hole symmetric 
ICPQMC simulation of Dagotto et al. 6] to systems 
without this symmetry, such as the tetrahedron, trun- 
cated tetrahedron, and Cgo molecule. Our simulations 
show that an accurate canonical partition function ratio 
can be obtained through this technique. Unfortunately, 
the fitting of these ratios to obtain accurate energy gaps 
for Ceo is impractical. Nevertheless, consistency between 
ED, PQMC and ICPQMC has been found. 
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